The packages are automatically loaded using pacman. The reported .html was last ran on the system: x86_64-apple-darwin17.0 and R version: R version 4.2.1 (2022-06-23)

1 Sim data

First, for the Stage 1 registered report some ICC distributions are simulated, aiming for a bimodial distribution to pick from and some specific categories (e.g., low, high, medium). Since rnorm will simulate beyond what is preferred, subseting < 1.00 and > -(.15).

1.1 simulate ICC categories

# three  diff ICC groups
# Elliott et al & Noble et al have largely reported a concentration around .3-.5. Here, simulate
# based on four levels, where min > -.15 and max is < 1.00
low_neg <- rnorm(10000, mean = 0, sd = 0.05)
low <- rnorm(10000, mean = 0.15, sd = 0.05) %>% subset(. < 1 & . > -.15)
mid <- rnorm(10000, mean = 0.45, sd = 0.1) %>% subset(. < 1 & . > -.15)
high <- rnorm(10000, mean = 0.7, sd = 0.2) %>% subset(. < 1 & . > -.15)
values <- c(low_neg, low, mid,high)

plotting histogram for all values

hist(values, breaks = 10, col = "gray", xlab = "Values", main = "Bimodal Distribution of ICCs")

1.2 defining & creating model permutations

Here, permutations are based on FWHM = five options, Motion = six options, Task Models = three options, Task Contrasts = four options, across three studies total 360 options that will be the “final” simulated dataset.

fwhm_ops = c(4,5,6,7,8) # five FWHM

motion_ops = c("None", "Regress_TransRot","Regress_Derivatives","Regress_1st8aCompCor",
               "Censor_mfd","Excl_mFD") # fd cutoff, => .9

taskmod_ops = c("Cue_CueDur","Cue_CueFixDur","Fix_FixDur")
taskcontrast_ops = c("BigWin_v_Neut","BigWin_v_Base","SmWin_v_Neut","SmWin_v_Base")

studies = c("MLS","ABCD","AHRB")

# permutations of the variables & setting up for loop
var_permutations <- as.matrix(expand.grid(
  study = studies, fwhm = fwhm_ops, motion = motion_ops, 
  modeltype = taskmod_ops, contrast = taskcontrast_ops))

Creating empty DF to store info in

df = data.frame(study = character(), fwhm = numeric(),
                motion = character(), contrast = character(), modeltype = character(),
                med_icc = numeric())

1.3 Specifying sim ICCs

Here, using an abstract method (i.e., several if/else rules) to specify from which category of ICCs to sample from for specific modeling options. Note, this doesnt result in a covariance within studies, so the within-subj variance is likely to be random thus zero.

set.seed(1000)
# simulate across study
for (row in 1:nrow(var_permutations)) {
  
  
  tmp_fwhm =    var_permutations[row,"fwhm"]
  tmp_motion =  var_permutations[row,"motion"]
  tmp_contrast =    var_permutations[row,"contrast"]
  tmp_model =    var_permutations[row,"modeltype"]
  tmp_study = var_permutations[row,"study"]
  
  df[row,"fwhm"] <- tmp_fwhm
  df[row,"motion"] <- tmp_motion
  df[row,"contrast"] <- tmp_contrast
  df[row,"modeltype"] <- tmp_model
  df[row,"study"] <- tmp_study

  if (tmp_fwhm %in% c(7,8)) {
    if (tmp_motion %in% c("Regress_TransRot", "Regress_Derivatives","Regress_1st8aCompCor")) {
      if (tmp_model %in% c("Cue_CueDur","Cue_CueFixDur")) {
        if (tmp_contrast=="BigWin_v_Neut") {
          df[row,"med_icc"] <- sample(x = high, size = 1,replace = FALSE)
        } else {
          df[row,"med_icc"] <- sample(x = mid, size = 1,replace = FALSE)
        }
      } else {
        df[row,"med_icc"] <- sample(x = values, size = 1,replace = FALSE)
      }
    } else {
      df[row,"med_icc"] <- sample(x = values, size = 1,replace = FALSE)
    }
  } else {
    if (tmp_motion %in% c("Regress_TransRot", "Regress_Derivatives","Regress_1st8aCompCor")) {
      if (tmp_model %in% c("Cue_CueDur","Cue_CueFixDur")) {
        if (tmp_contrast=="BigWin_v_Neut") {
          df[row,"med_icc"] <- sample(x = mid, size = 1,replace = FALSE)
        } else {
          df[row,"med_icc"] <- sample(x = low, size = 1,replace = FALSE)
        }
      } else {
        df[row,"med_icc"] <- sample(x = low, size = 1,replace = FALSE)
      }
    } else {
      df[row,"med_icc"] <- sample(x = low_neg, size = 1,replace = FALSE)
    }
  }
}

releving and specifying factors so things are interpretable in the heiarchical models.

# level/relabel
df$fwhm_r <- as.numeric(df$fwhm)
df$motion = relevel(as.factor(df$motion), ref = "None", 
                     levels = c("None","Censor_mfd","Excl_mFD",
                                "Regress_1st8aCompCor","Regress_Derivatives",
                                "Regress_TransRot"))
df$modeltype = relevel(as.factor(df$modeltype), ref = "Cue_CueDur", 
                     levels = c("Cue_CueDur","Cue_CuFixDur","Fix_FixDur"))
df$contrast = relevel(as.factor(df$contrast), ref = "BigWin_v_Base", 
                     levels = c("BigWin_v_Base","BigWin_v_Neut","SmWin_v_Base","SmWin_v_Neut"))

2 Results: ICCs across analytic approaches

Below, running the steps to summarize the different ICCs from the model combinations 360 for each study

2.1 plot distribution

Plotting overall and for each of [four] categories

icc_dist = df %>% 
  ggplot(aes(x = med_icc)) + 
  geom_histogram(bins = 50, fill = "white", colour = "black") +
  ggtitle("Overall Distribution") +
  theme_minimal()


fwhm_dist = df %>% 
  ggplot(aes(x = med_icc)) +
  geom_histogram(bins = 50, fill = "white", colour = "black") +
  facet_wrap(~fwhm) +
  ggtitle("Distribution by FWHM category") +
  theme_minimal()


motion_dist = df %>% 
  ggplot(aes(x = med_icc)) +
  geom_histogram(bins = 50, fill = "white", colour = "black") +
  facet_wrap(~motion) +
  ggtitle("Distribution by Motion category") +
  theme_minimal()

contrast_dist = df %>% 
  ggplot(aes(x = med_icc)) +
  geom_histogram(bins = 50, fill = "white", colour = "black") +
  facet_wrap(~contrast) +
  ggtitle("Distribution by Contrast category") +
  theme_minimal()

modeltype_dist = df %>% 
  ggplot(aes(x = med_icc)) +
  geom_histogram(bins = 50, fill = "white", colour = "black") +
  facet_wrap(~modeltype) +
  ggtitle("Distribution by Model Type category") +
  theme_minimal()

icc_dist

fwhm_dist

motion_dist

contrast_dist

modeltype_dist

Creating rainclouds via ggrain

fwhm_rg = df %>% ggplot(aes(x = fwhm, y = med_icc, fill = fwhm, color = fwhm)) +
  geom_rain(alpha = .5, rain.side = 'l',
            boxplot.args = list(color = "black", outlier.shape = NA),
            boxplot.args.pos = list(
              position = ggpp::position_dodgenudge(x = .1), width = 0.1
            )) +
  theme_classic() +
  ggtitle("Distribution by FWHN category") +
  scale_fill_brewer(palette = 'Dark2') +
  scale_color_brewer(palette = 'Dark2') +
  guides(fill = 'none', color = 'none') +
  theme(text = element_text(family = "Times New Roman"),
        axis.text = element_text(size = 10),
        axis.title = element_text(size = 10),
        legend.text = element_text(size = 10),
        legend.title = element_text(size = 10),
        plot.title = element_text(size = 16))


motion_rg = df %>% ggplot(aes(x = motion, y = med_icc, fill = motion, color = motion)) +
  geom_rain(alpha = .5, rain.side = 'l',
            boxplot.args = list(color = "black", outlier.shape = NA),
            boxplot.args.pos = list(
              position = ggpp::position_dodgenudge(x = .1), width = 0.1
            )) +
  theme_classic() +
  ggtitle("Distribution by Motion category") +
  scale_fill_brewer(palette = 'Dark2') +
  scale_color_brewer(palette = 'Dark2') +
  guides(fill = 'none', color = 'none') +
  theme(text = element_text(family = "Times New Roman"),
        axis.text = element_text(size = 10),
        axis.title = element_text(size = 10),
        legend.text = element_text(size = 10),
        legend.title = element_text(size = 10),
        plot.title = element_text(size = 16))

modeltype_rg = df %>% ggplot(aes(x = modeltype, y = med_icc, fill = modeltype, color = modeltype)) +
  geom_rain(alpha = .5, rain.side = 'l',
            boxplot.args = list(color = "black", outlier.shape = NA),
            boxplot.args.pos = list(
              position = ggpp::position_dodgenudge(x = .1), width = 0.1
            )) +
  theme_classic() +
  ggtitle("Distribution by Model Type category") +
  scale_fill_brewer(palette = 'Dark2') +
  scale_color_brewer(palette = 'Dark2') +
  guides(fill = 'none', color = 'none') +
  theme(text = element_text(family = "Times New Roman"),
        axis.text = element_text(size = 10),
        axis.title = element_text(size = 10),
        legend.text = element_text(size = 10),
        legend.title = element_text(size = 10),
        plot.title = element_text(size = 16))

contrast_rg = df %>% ggplot(aes(x = contrast, y = med_icc, fill = contrast, color = contrast)) +
  geom_rain(alpha = .5, rain.side = 'l',
            boxplot.args = list(color = "black", outlier.shape = NA),
            boxplot.args.pos = list(
              position = ggpp::position_dodgenudge(x = .1), width = 0.1
            )) +
  theme_classic() +
  ggtitle("Distribution by Contrast category") +
  scale_fill_brewer(palette = 'Dark2') +
  scale_color_brewer(palette = 'Dark2') +
  guides(fill = 'none', color = 'none') +
  theme(text = element_text(family = "Times New Roman"),
        axis.text = element_text(size = 10),
        axis.title = element_text(size = 10),
        legend.text = element_text(size = 10),
        legend.title = element_text(size = 10),
        plot.title = element_text(size = 16))

fwhm_rg

motion_rg

modeltype_rg

contrast_rg

2.2 summaries of model types

summary_icc = df %>% 
  group_by(fwhm, motion, modeltype, contrast) %>% 
  summarise(across(("med_icc"), list(min = min, max = max, mean = mean, median = median))) %>% 
  as.data.frame()

## `summarise()` has grouped output by 'fwhm', 'motion', 'modeltype'. You can
## override using the `.groups` argument.

summary_icc %>% 
  arrange(desc(med_icc_min)) %>% 
  slice(1:10)

##    fwhm               motion     modeltype      contrast med_icc_min
## 1     7 Regress_1st8aCompCor Cue_CueFixDur BigWin_v_Neut   0.7562752
## 2     7     Regress_TransRot    Cue_CueDur BigWin_v_Neut   0.7342187
## 3     7 Regress_1st8aCompCor    Cue_CueDur BigWin_v_Neut   0.6781959
## 4     8 Regress_1st8aCompCor    Cue_CueDur BigWin_v_Neut   0.6599825
## 5     8     Regress_TransRot Cue_CueFixDur BigWin_v_Neut   0.6202459
## 6     8     Regress_TransRot    Cue_CueDur BigWin_v_Neut   0.6006907
## 7     7  Regress_Derivatives    Cue_CueDur BigWin_v_Neut   0.5129165
## 8     7  Regress_Derivatives    Cue_CueDur  SmWin_v_Base   0.5122227
## 9     8     Regress_TransRot    Cue_CueDur  SmWin_v_Base   0.4910063
## 10    7           Censor_mfd Cue_CueFixDur BigWin_v_Base   0.4752819
##    med_icc_max med_icc_mean med_icc_median
## 1    0.9152130    0.8166406      0.7784335
## 2    0.8155691    0.7797988      0.7896087
## 3    0.8508815    0.7753146      0.7968663
## 4    0.8927606    0.7572947      0.7191411
## 5    0.8024149    0.7262381      0.7560534
## 6    0.7947373    0.6947145      0.6887155
## 7    0.6764460    0.5986463      0.6065764
## 8    0.6208878    0.5716596      0.5818682
## 9    0.5377776    0.5122792      0.5080539
## 10   0.7166981    0.6078405      0.6315416

summary_icc %>% 
  arrange(desc(med_icc_max)) %>% 
  slice(1:10)

##    fwhm               motion     modeltype      contrast med_icc_min
## 1     7           Censor_mfd    Fix_FixDur  SmWin_v_Base  0.02215935
## 2     7           Censor_mfd    Fix_FixDur BigWin_v_Base  0.02659101
## 3     8             Excl_mFD Cue_CueFixDur  SmWin_v_Neut  0.10338313
## 4     7             Excl_mFD    Cue_CueDur BigWin_v_Base  0.08332403
## 5     7     Regress_TransRot    Fix_FixDur  SmWin_v_Neut  0.41721161
## 6     7 Regress_1st8aCompCor Cue_CueFixDur BigWin_v_Neut  0.75627515
## 7     8             Excl_mFD    Cue_CueDur BigWin_v_Neut  0.12289872
## 8     8                 None    Cue_CueDur  SmWin_v_Neut -0.04767462
## 9     8 Regress_1st8aCompCor    Cue_CueDur BigWin_v_Neut  0.65998255
## 10    7                 None    Fix_FixDur  SmWin_v_Neut -0.05309287
##    med_icc_max med_icc_mean med_icc_median
## 1    0.9947495    0.6162410      0.8318142
## 2    0.9640912    0.4915768      0.4840481
## 3    0.9435615    0.3880643      0.1172483
## 4    0.9398068    0.4413431      0.3008984
## 5    0.9364607    0.6109669      0.4792284
## 6    0.9152130    0.8166406      0.7784335
## 7    0.9014496    0.5025672      0.4833531
## 8    0.8950570    0.4689903      0.5595883
## 9    0.8927606    0.7572947      0.7191411
## 10   0.8536519    0.3142251      0.1421162

summary_icc %>% 
  arrange(desc(med_icc_mean)) %>% 
  slice(1:10)

##    fwhm               motion     modeltype      contrast med_icc_min
## 1     7 Regress_1st8aCompCor Cue_CueFixDur BigWin_v_Neut  0.75627515
## 2     7     Regress_TransRot    Cue_CueDur BigWin_v_Neut  0.73421866
## 3     7 Regress_1st8aCompCor    Cue_CueDur BigWin_v_Neut  0.67819590
## 4     8 Regress_1st8aCompCor    Cue_CueDur BigWin_v_Neut  0.65998255
## 5     8     Regress_TransRot Cue_CueFixDur BigWin_v_Neut  0.62024593
## 6     8     Regress_TransRot    Cue_CueDur BigWin_v_Neut  0.60069068
## 7     8 Regress_1st8aCompCor Cue_CueFixDur BigWin_v_Neut  0.42787099
## 8     8  Regress_Derivatives Cue_CueFixDur BigWin_v_Neut  0.44973463
## 9     7  Regress_Derivatives Cue_CueFixDur BigWin_v_Neut  0.44215137
## 10    7           Censor_mfd    Fix_FixDur  SmWin_v_Base  0.02215935
##    med_icc_max med_icc_mean med_icc_median
## 1    0.9152130    0.8166406      0.7784335
## 2    0.8155691    0.7797988      0.7896087
## 3    0.8508815    0.7753146      0.7968663
## 4    0.8927606    0.7572947      0.7191411
## 5    0.8024149    0.7262381      0.7560534
## 6    0.7947373    0.6947145      0.6887155
## 7    0.8232779    0.6800934      0.7891313
## 8    0.8360832    0.6592579      0.6919558
## 9    0.7777469    0.6316792      0.6751393
## 10   0.9947495    0.6162410      0.8318142

2.3 Fitting HLM

Below using lmer to fit HLM model. One model has random intercept (study), second model has random intercept (study) and random slope (FWHM). Reporting using tab_model

Level 1 model

\(MedianICC_{ij} = \beta_{0_j} + \beta_{1_j} \times \text{fwhm}_{ij} + \beta_{2_j} \times \text{motion}_{ij} + \beta_{3_j} \times \text{modeltype}_{ij} + \beta_{4_j} \times \text{contrast}_{ij} + \epsilon_{ij}\)

Level 2 model

\(\beta_{0_j} = \gamma_{00} + u_{0_j}\)

\(\beta_{1_j} = \gamma_{10} + u_{1_j}\)

\(\beta_{2_j} = \gamma_{20} + u_{2_j}\)

\(\beta_{3_j} = \gamma_{30} + u_{3_j}\)

\(\beta_{4_j} = \gamma_{40} + u_{4_j}\)

`Notes:

med_icc_ij is the outcome for the ith model in the jth sample β0_j is the intercept for the jth sample fwhm_r_ij, motion_ij, modeltype_ij, and contrast_ij are the predictor variables for the ith model in the jth sample β1_j, β2_j, β3_j, and β4_j: corresponding coefficients for the jth sample ε_ij is the within-sample error term

γ00, γ10, γ20, γ30, and γ40 are the fixed effects intercept and slopes for the predictors; u0_j, u1_j, u2_j, u3_j, and u4_j are the sample-level residuals (random effects) for the intercept and the four slopes.`

model = lmer(med_icc ~ fwhm_r + motion + modeltype + contrast + (1 | study), data = df)
model_rndslp = lmer(med_icc ~ fwhm_r + motion + modeltype + contrast + (1+fwhm_r| study), data = df)

## boundary (singular) fit: see help('isSingular')

tab_model(model, model_rndslp)

	med_icc			med_icc
Predictors	Estimates	CI	p	Estimates	CI	p
(Intercept)	-0.35	-0.41 – -0.28	<0.001	-0.35	-0.42 – -0.27	<0.001
fwhm r	0.08	0.07 – 0.09	<0.001	0.08	0.07 – 0.09	<0.001
motion [Censor_mfd]	-0.02	-0.06 – 0.02	0.382	-0.02	-0.06 – 0.02	0.381
motion [Excl_mFD]	0.01	-0.03 – 0.05	0.667	0.01	-0.03 – 0.05	0.667
motion [Regress_1st8aCompCor]	0.16	0.12 – 0.20	<0.001	0.16	0.12 – 0.20	<0.001
motion [Regress_Derivatives]	0.15	0.11 – 0.18	<0.001	0.15	0.11 – 0.18	<0.001
motion [Regress_TransRot]	0.17	0.13 – 0.21	<0.001	0.17	0.13 – 0.21	<0.001
modeltype [Cue_CueFixDur]	-0.01	-0.04 – 0.02	0.432	-0.01	-0.04 – 0.02	0.431
modeltype [Fix_FixDur]	-0.08	-0.10 – -0.05	<0.001	-0.08	-0.10 – -0.05	<0.001
contrast [BigWin_v_Neut]	0.08	0.05 – 0.12	<0.001	0.08	0.05 – 0.12	<0.001
contrast [SmWin_v_Base]	0.00	-0.03 – 0.03	0.991	0.00	-0.03 – 0.03	0.991
contrast [SmWin_v_Neut]	-0.01	-0.04 – 0.03	0.721	-0.01	-0.04 – 0.03	0.720
Random Effects
σ²	0.04			0.04
τ₀₀	0.00 _study			0.00 _study
τ₁₁				0.00 _{study.fwhm_r}
ρ₀₁				-1.00 _study
ICC	0.00
N	3 _study			3 _study
Observations	1080			1080
Marginal R² / Conditional R²	0.383 / 0.383			0.384 / NA

Using emmeans to control type I error rate of controls via Tukey’s Honest Significant Test

# Model comparisons
mod1_lmer_obs = as.numeric(summary(model)[3]$devcomp$dims[1])
emm_options(lmerTest.limit = mod1_lmer_obs)

# contrast contrast
mod1_thsd = emmeans(object = model, specs = pairwise ~ contrast, pbkrtest.limit = mod1_lmer_obs)
mod1_thsd$contrasts %>% summary(infer = TRUE) %>% 
  kbl(digits = 4, caption = "Contrast Tukey HSD: Contrast") %>%  kable_styling()

Contrast Tukey HSD: Contrast
contrast	estimate	SE	df	lower.CL	upper.CL	t.ratio	p.value
BigWin_v_Base - BigWin_v_Neut	-0.0846	0.0163	1066	-0.1265	-0.0426	-5.1852	0.0000
BigWin_v_Base - SmWin_v_Base	-0.0002	0.0163	1066	-0.0421	0.0418	-0.0110	1.0000
BigWin_v_Base - SmWin_v_Neut	0.0058	0.0163	1066	-0.0361	0.0478	0.3578	0.9843
BigWin_v_Neut - SmWin_v_Base	0.0844	0.0163	1066	0.0424	0.1264	5.1742	0.0000
BigWin_v_Neut - SmWin_v_Neut	0.0904	0.0163	1066	0.0484	0.1324	5.5430	0.0000
SmWin_v_Base - SmWin_v_Neut	0.0060	0.0163	1066	-0.0360	0.0480	0.3688	0.9829

# contrast motion
mod1_thsd = emmeans(object = model, specs = pairwise ~ motion, pbkrtest.limit = mod1_lmer_obs)
mod1_thsd$contrasts %>% summary(infer = TRUE) %>% 
  kbl(digits = 4, caption = "Contrast Tukey HSD: Motion Corr") %>%  kable_styling()

Contrast Tukey HSD: Motion Corr
contrast	estimate	SE	df	lower.CL	upper.CL	t.ratio	p.value
None - Censor_mfd	0.0175	0.02	1066	-0.0396	0.0745	0.8738	0.9527
None - Excl_mFD	-0.0086	0.02	1066	-0.0656	0.0484	-0.4297	0.9981
None - Regress_1st8aCompCor	-0.1605	0.02	1066	-0.2175	-0.1035	-8.0352	0.0000
None - Regress_Derivatives	-0.1455	0.02	1066	-0.2026	-0.0885	-7.2862	0.0000
None - Regress_TransRot	-0.1666	0.02	1066	-0.2236	-0.1095	-8.3391	0.0000
Censor_mfd - Excl_mFD	-0.0260	0.02	1066	-0.0831	0.0310	-1.3035	0.7832
Censor_mfd - Regress_1st8aCompCor	-0.1780	0.02	1066	-0.2350	-0.1209	-8.9089	0.0000
Censor_mfd - Regress_Derivatives	-0.1630	0.02	1066	-0.2200	-0.1060	-8.1599	0.0000
Censor_mfd - Regress_TransRot	-0.1840	0.02	1066	-0.2410	-0.1270	-9.2128	0.0000
Excl_mFD - Regress_1st8aCompCor	-0.1519	0.02	1066	-0.2089	-0.0949	-7.6054	0.0000
Excl_mFD - Regress_Derivatives	-0.1370	0.02	1066	-0.1940	-0.0799	-6.8564	0.0000
Excl_mFD - Regress_TransRot	-0.1580	0.02	1066	-0.2150	-0.1010	-7.9093	0.0000
Regress_1st8aCompCor - Regress_Derivatives	0.0150	0.02	1066	-0.0421	0.0720	0.7490	0.9756
Regress_1st8aCompCor - Regress_TransRot	-0.0061	0.02	1066	-0.0631	0.0510	-0.3039	0.9997
Regress_Derivatives - Regress_TransRot	-0.0210	0.02	1066	-0.0781	0.0360	-1.0529	0.8997

# contrast modeltype
mod1_thsd = emmeans(object = model, specs = pairwise ~ modeltype, pbkrtest.limit = mod1_lmer_obs)
mod1_thsd$contrasts %>% summary(infer = TRUE) %>% 
  kbl(digits = 4, caption = "Contrast Tukey HSD: Model Type ") %>%  kable_styling()

Contrast Tukey HSD: Model Type
contrast	estimate	SE	df	lower.CL	upper.CL	t.ratio	p.value
Cue_CueDur - Cue_CueFixDur	0.0111	0.0141	1066	-0.0220	0.0442	0.7858	0.7118
Cue_CueDur - Fix_FixDur	0.0764	0.0141	1066	0.0433	0.1096	5.4101	0.0000
Cue_CueFixDur - Fix_FixDur	0.0653	0.0141	1066	0.0322	0.0985	4.6242	0.0000

2.4 Specification Curve

Creating data in a format that is compatible with specr. Needs: estimate (i.e., ICC), std.error, conf.high, conf.low.

# create specr plot for med_icc averages for 
# first, combine independent model vars into string to create average for each model type
df$model_type <- paste(df$fwhm,df$motion,df$contrast,df$modeltype,sep = "-")

# calculate the avg estimate of ICC across study, standard error and +/- 95% confidence interval. In complete version
df_summ <- df %>% 
  group_by(model_type) %>% 
  summarise(estimate = mean(med_icc), std.error = sd(med_icc)/sqrt(length(med_icc))) %>% 
  mutate(conf.low = estimate - 1.96 * std.error, conf.high = estimate + 1.96 * std.error) %>% 
  separate(col = model_type, into = c("fwhm","motion","contrast","modeltype"), sep = "-", remove = FALSE)

2.4.1 All models

creating combined panel 1 and panel 2 for all model permutations first.

plot_a = plot_curve(df = df_summ, ci = TRUE, desc = FALSE, legend = FALSE, null = 0)

plot_b <- plot_choices(df = df_summ, choices = c("fwhm", "motion","contrast","modeltype"), desc = F, null = 0) +
  labs(y = "Variables", x = "Ordered Specification Curve \n ICC coefficient")

cowplot::plot_grid(plot_a, plot_b, ncol = 1, align = "v", axis = 'tblr',
                   labels = c('A', 'B'), rel_heights = c(1, 2),
                   label_fontfamily = "Times", label_size = 12)

2.4.2 subset models

Creating model that is subset to visualize reliability for the top (>75th) and bottom (< 25th) quartile

# get 75th/25th qunatiles
top_75q = as.numeric(quantile(df_summ$estimate, .75))
bot_25q = as.numeric(quantile(df_summ$estimate, .25))
df_summ_subset = df_summ %>% 
  filter(estimate < bot_25q | estimate > top_75q)

# subset plots
plot_a_sub = plot_curve(df = df_summ_subset, ci = TRUE, desc = FALSE, legend = FALSE, null = 0)

plot_b_sub <- plot_choices(df = df_summ_subset, choices = c("fwhm", "motion","contrast","modeltype"), desc = F, null = 0) +
  labs(y = "Variables", x = "Ordered Specification Curve \n ICC coefficient")

cowplot::plot_grid(plot_a_sub, plot_b_sub, ncol = 1, align = "v", axis = 'tblr',
                   labels = c('A', 'B'), rel_heights = c(1, 2),
                   label_fontfamily = "Times", label_size = 12)

3 Design efficiency Examples

The ability to estimate the BOLD signal in the timeseries, efficiency plays a significant role. See issues discussed by series of videos from Jeanette Mumford. The MID task, has historically different ways of defining the “Anticipation” phase. Previously, it started by modeling the fixation phase/duration and later some have used the cue/duration. However, the issues have been rarely cleared up in the literature as it is challenging to decipher from the information provided in papers what is defined as the anticipation phase and duration. Hence, up to three models can be created in the design matrix.

Below, using the efficiency calculation tutorial from Dr Mumford to estimate the effieiency of the three different models that can be estimated in this version of the MID task.

3.1 Pulling Behavior Data

First, the behavior data is pulled from the events files for each subject and each run. Some subjects are excluded because they lack hit/miss trials which causes an error in the step using neuRosim.

# location for subjects run behavior data
# update path as needed
dir = "/Users/michaeldemidenko/Desktop/Academia/UM/2_AHRB/Projects/BIDS/Beh/MID_W1/"

# getting list of file names
beh_files = list.files(dir)

# removing subjects that dont have some Hit/Miss for outcome, as there are zero values for complete calculation
beh_files_ct = beh_files[-c(12,19,34,39,46,120,140)]
#   sub-07_ses-1_task-mid_run-02_events.tsv
#   sub-100_ses-1_task-mid_run-01_events.tsv
#   sub-107_ses-1_task-mid_run-02_events.tsv
#   sub-12_ses-1_task-mid_run-01_events.tsv
#   sub-15_ses-1_task-mid_run-02_events.tsv
#   sub-52_ses-1_task-mid_run-02_events.tsv
#   sub-62_ses-1_task-mid_run-02_events.tsv

3.2 Estimating efficiency from Events

Using for loop to create the the data based on the behavioral onsets/durations for efficiency results. Three different models are fit:

design_CueMod: Cue Onset + Cue Duration
design_AntMod: Cue Onset + (Cue Duration + Fixation Duration)
design_FixMod: Fixation Onset + Fixation Duration

While the focus is on the anticipatory phase, e.g. modeling 5 cue times + associated duration, the feedback phase is modeled, too. In other words, the total model includes 15 items.

# Null df to add to
df = NULL

# creating forloop using Jeanette Mumfords recommended code from Neurohack - https://github.com/jmumford/nhwEfficiency
# tailored to work with behavior data collected for MID task data

#set TR, vols, effect size (amplitude) and run and mri length
tr = .8 
vols = 407
effect_size = 1
run_len = vols*tr
mri_len = seq(1,407,1)
  
for (f in beh_files_ct) {
  
  # Combining directory and file path to open .csv
  beh_dat = read.csv(paste(dir,f,sep=""),sep = '\t')
  
  # subj/run labels to use in later df creation
  str_vals =data.frame(str_split(f,'_'))
  sub = str_vals[1,1]
  run = str_vals[4,1]
  
  # convert one of the outcome items foor Neutral condition so it is differently labeled for Hit/Miss
  beh_dat <- beh_dat %>% 
    mutate(TRIAL_RESULT = if_else(TRIAL_RESULT!="No money at stake!",TRIAL_RESULT,
                                  if_else(PROBE_HIT==1,TRIAL_RESULT,"No Money At Stake, MISS")))
  
  #
  
  ### Anticipation
  
  # Onsets
  BW_ons_1 <-   beh_dat$CUE_ONSET[beh_dat$TRIAL_TYPE=="LargeGain"]
  SW_ons_1 <-   beh_dat$CUE_ONSET[beh_dat$TRIAL_TYPE=="SmallGain"]
  Neut_ons_1 <- beh_dat$CUE_ONSET[beh_dat$TRIAL_TYPE=="NoMoneyStake"]
  BL_ons_1 <-   beh_dat$CUE_ONSET[beh_dat$TRIAL_TYPE=="LargeLoss"]
  SL_ons_1 <-   beh_dat$CUE_ONSET[beh_dat$TRIAL_TYPE=="SmallLoss"]
  
  BW_ons_2 <-   beh_dat$FIXATION_ONSET[beh_dat$TRIAL_TYPE=="LargeGain"]
  SW_ons_2 <-   beh_dat$FIXATION_ONSET[beh_dat$TRIAL_TYPE=="SmallGain"]
  Neut_ons_2 <- beh_dat$FIXATION_ONSET[beh_dat$TRIAL_TYPE=="NoMoneyStake"]
  BL_ons_2 <-   beh_dat$FIXATION_ONSET[beh_dat$TRIAL_TYPE=="LargeLoss"]
  SL_ons_2 <-   beh_dat$FIXATION_ONSET[beh_dat$TRIAL_TYPE=="SmallLoss"]
  
  # durations
  BW_dur1 <-   beh_dat$CUE_DURATION[beh_dat$TRIAL_TYPE=="LargeGain"]   
  SW_dur1 <-   beh_dat$CUE_DURATION[beh_dat$TRIAL_TYPE=="SmallGain"]   
  Neut_dur1 <- beh_dat$CUE_DURATION[beh_dat$TRIAL_TYPE=="NoMoneyStake"]
  BL_dur1 <-   beh_dat$CUE_DURATION[beh_dat$TRIAL_TYPE=="LargeLoss"]   
  SL_dur1 <-   beh_dat$CUE_DURATION[beh_dat$TRIAL_TYPE=="SmallLoss"]   
  
  BW_dur2 <-   beh_dat$FIXATION_DURATION[beh_dat$TRIAL_TYPE=="LargeGain"]   
  SW_dur2 <-   beh_dat$FIXATION_DURATION[beh_dat$TRIAL_TYPE=="SmallGain"]   
  Neut_dur2 <- beh_dat$FIXATION_DURATION[beh_dat$TRIAL_TYPE=="NoMoneyStake"]
  BL_dur2 <-   beh_dat$FIXATION_DURATION[beh_dat$TRIAL_TYPE=="LargeLoss"]   
  SL_dur2 <-   beh_dat$FIXATION_DURATION[beh_dat$TRIAL_TYPE=="SmallLoss"]   
  
  BW_dur12 <-   beh_dat$CUE_DURATION[beh_dat$TRIAL_TYPE=="LargeGain"]   + beh_dat$FIXATION_DURATION[beh_dat$TRIAL_TYPE=="LargeGain"]   
  SW_dur12 <-   beh_dat$CUE_DURATION[beh_dat$TRIAL_TYPE=="SmallGain"]   + beh_dat$FIXATION_DURATION[beh_dat$TRIAL_TYPE=="SmallGain"]   
  Neut_dur12 <- beh_dat$CUE_DURATION[beh_dat$TRIAL_TYPE=="NoMoneyStake"]+ beh_dat$FIXATION_DURATION[beh_dat$TRIAL_TYPE=="NoMoneyStake"]
  BL_dur12 <-   beh_dat$CUE_DURATION[beh_dat$TRIAL_TYPE=="LargeLoss"]   + beh_dat$FIXATION_DURATION[beh_dat$TRIAL_TYPE=="LargeLoss"]   
  SL_dur12 <-   beh_dat$CUE_DURATION[beh_dat$TRIAL_TYPE=="SmallLoss"]   + beh_dat$FIXATION_DURATION[beh_dat$TRIAL_TYPE=="SmallLoss"]   
  
  
  ### Feedback
  # Onsets
  BW_Hit_ons <-   beh_dat$FEEDBACK_ONSET[beh_dat$TRIAL_RESULT=="You earn $5!"]
  SW_Hit_ons <-   beh_dat$FEEDBACK_ONSET[beh_dat$TRIAL_RESULT=="You earn $0.20!"]
  Neut_Hit_ons <- beh_dat$FEEDBACK_ONSET[beh_dat$TRIAL_RESULT=="No money at stake!"]
  BL_Hit_ons <-   beh_dat$FEEDBACK_ONSET[beh_dat$TRIAL_RESULT=="You keep $5!"]
  SL_Hit_ons <-   beh_dat$FEEDBACK_ONSET[beh_dat$TRIAL_RESULT=="You keep $0.20!"]
  
  BW_Miss_ons <-   beh_dat$FEEDBACK_ONSET[beh_dat$TRIAL_RESULT=="You did not earn $5!"]
  SW_Miss_ons <-   beh_dat$FEEDBACK_ONSET[beh_dat$TRIAL_RESULT=="You did not earn $0.20!"]
  Neut_Miss_ons <- beh_dat$FEEDBACK_ONSET[beh_dat$TRIAL_RESULT=="No Money At Stake, MISS"]
  BL_Miss_ons <-   beh_dat$FEEDBACK_ONSET[beh_dat$TRIAL_RESULT=="You lose $5!"]
  SL_Miss_ons <-   beh_dat$FEEDBACK_ONSET[beh_dat$TRIAL_RESULT=="You lose $0.20!"]
  
  # durations
  BW_Hit_dur <-   beh_dat$FEEDBACK_DURATION[beh_dat$TRIAL_RESULT=="You earn $5!"]
  SW_Hit_dur <-   beh_dat$FEEDBACK_DURATION[beh_dat$TRIAL_RESULT=="You earn $0.20!"]
  Neut_Hit_dur <- beh_dat$FEEDBACK_DURATION[beh_dat$TRIAL_RESULT=="No money at stake!"]
  BL_Hit_dur <-   beh_dat$FEEDBACK_DURATION[beh_dat$TRIAL_RESULT=="You keep $5!"]
  SL_Hit_dur <-   beh_dat$FEEDBACK_DURATION[beh_dat$TRIAL_RESULT=="You keep $0.20!"]
  
  BW_Miss_dur <-   beh_dat$FEEDBACK_DURATION[beh_dat$TRIAL_RESULT=="You did not earn $5!"]
  SW_Miss_dur <-   beh_dat$FEEDBACK_DURATION[beh_dat$TRIAL_RESULT=="You did not earn $0.20!"]
  Neut_Miss_dur <- beh_dat$FEEDBACK_DURATION[beh_dat$TRIAL_RESULT=="No Money At Stake, MISS"]
  BL_Miss_dur <-   beh_dat$FEEDBACK_DURATION[beh_dat$TRIAL_RESULT=="You lose $5!"]
  SL_Miss_dur <-   beh_dat$FEEDBACK_DURATION[beh_dat$TRIAL_RESULT=="You lose $0.20!"]
  
  
  ## model 1 (cue only), creating a list of values for Onset and Duration
  CueMod_ons = list(BW_ons_1, SW_ons_1,Neut_ons_1,BL_ons_1,SL_ons_1,
                  BW_Hit_ons, SW_Hit_ons, Neut_Hit_ons, BL_Hit_ons,SL_Hit_ons,
                  BW_Miss_ons,SW_Miss_ons, Neut_Miss_ons, BL_Miss_ons,SL_Miss_ons
                  )
  CueMod_dur = list(BW_dur1,SW_dur1,Neut_dur1,BL_dur1,SL_dur1,
                  BW_Hit_dur, SW_Hit_dur, Neut_Hit_dur, BL_Hit_dur,SL_Hit_dur,
                  BW_Miss_dur,SW_Miss_dur, Neut_Miss_dur, BL_Miss_dur,SL_Miss_dur
                  )
  
  # Using specifydesign to create the HRF Cue Model based on inputs:
  # onsets of the specified model
  # durations of the specificated model
  # the list of effect size (or amplitudes) used for each stimulus in model
  # Tr and total time to create the length of HRF using the convolved double-gamma
  # condition names is a list of labels for the Input onset labels
  Cue_Mod <- specifydesign(onsets = CueMod_ons,
                      durations = CueMod_dur, 
                      effectsize =list(1,1,1,1,1,
                                       1,1,1,1,1,
                                       1,1,1,1,1),
                      TR = tr, totaltime = run_len, conv = "double-gamma",
                      cond.names = list("BW","SW","Neut","BL","SL",
                                        "BWhit","SWhit","NEUThit","BLhit","SLhit",
                                        "BWmss","SWmss","NEUTmss","BLmss","SLmss")
                      )
  
  
  # Model 2 (fixation only)
  ## model 1
  FixMod_ons = list(BW_ons_2, SW_ons_2,Neut_ons_2,BL_ons_2,SL_ons_2,
                    BW_Hit_ons, SW_Hit_ons, Neut_Hit_ons, BL_Hit_ons,SL_Hit_ons,
                    BW_Miss_ons,SW_Miss_ons, Neut_Miss_ons, BL_Miss_ons,SL_Miss_ons
  )
  FixMod_dur = list(BW_dur2,SW_dur2,Neut_dur2,BL_dur2,SL_dur2,
                    BW_Hit_dur, SW_Hit_dur, Neut_Hit_dur, BL_Hit_dur,SL_Hit_dur,
                    BW_Miss_dur,SW_Miss_dur, Neut_Miss_dur, BL_Miss_dur,SL_Miss_dur
  )
  
  
  Fix_Mod <- specifydesign(onsets = FixMod_ons,
                           durations = FixMod_dur, 
                           effectsize =list(1,1,1,1,1,
                                            1,1,1,1,1,
                                            1,1,1,1,1),
                           TR = tr, totaltime = run_len, conv = "double-gamma",
                           cond.names = list("BW","SW","Neut","BL","SL",
                                             "BWhit","SWhit","NEUThit","BLhit","SLhit",
                                             "BWmss","SWmss","NEUTmss","BLmss","SLmss")
  )
  
  #Anticipation model cue onset & cue + fix dur
  ## model 1
  AntMod_ons = list(BW_ons_1, SW_ons_1,Neut_ons_1,BL_ons_1,SL_ons_1,
                    BW_Hit_ons, SW_Hit_ons, Neut_Hit_ons, BL_Hit_ons,SL_Hit_ons,
                    BW_Miss_ons,SW_Miss_ons, Neut_Miss_ons, BL_Miss_ons,SL_Miss_ons
  )
  AntMod_dur = list(BW_dur12,SW_dur12,Neut_dur12,BL_dur12,SL_dur12,
                    BW_Hit_dur, SW_Hit_dur, Neut_Hit_dur, BL_Hit_dur,SL_Hit_dur,
                    BW_Miss_dur,SW_Miss_dur, Neut_Miss_dur, BL_Miss_dur,SL_Miss_dur
  )
  
  
  Ant_Mod <- specifydesign(onsets = AntMod_ons,
                           durations = AntMod_dur, 
                           effectsize =list(1,1,1,1,1,
                                            1,1,1,1,1,
                                            1,1,1,1,1),
                           TR = tr, totaltime = run_len, conv = "double-gamma",
                           cond.names = list("BW","SW","Neut","BL","SL",
                                             "BWhit","SWhit","NEUThit","BLhit","SLhit",
                                             "BWmss","SWmss","NEUTmss","BLmss","SLmss")
  )
  
  design_CueMod <- data.frame("TR" = mri_len, Cue_Mod)

    # create design matrix + test contrast
  CueMod_mat = cbind("Inter" = rep(1,length(design_CueMod$BW)),
                     design_CueMod[2:16])
  # Plot fixation model and create design matrix
  design_FixMod <- data.frame("TR" = mri_len, Fix_Mod)
  
  # create design matrix 
  FixMod_mat = cbind("Inter" = rep(1,length(design_FixMod$BW)),
                     design_FixMod[2:16])
  
  # Plot anticipation model and create design matrix
  design_AntMod <- data.frame("TR" = mri_len, Ant_Mod)
  

  
  # create design matrix 
  AntMod_mat = cbind("Inter" = rep(1,length(design_AntMod$BW)),
                     design_AntMod[2:16])
  
  # setting up contrasts to estimate efficiency
  
  LGain = c(0,
          1,0,-1,0,0, # large gain versus neutral
          0,0,0,0,0,0,0,0,0,0
          )
  
  LGain_CueMod = 1/(t(LGain)%*%solve(t(CueMod_mat)%*%as.matrix(CueMod_mat))%*%LGain)
  LGain_AntMod = 1/(t(LGain)%*%solve(t(AntMod_mat)%*%as.matrix(AntMod_mat))%*%LGain)
  LGain_FixMod = 1/(t(LGain)%*%solve(t(FixMod_mat)%*%as.matrix(FixMod_mat))%*%LGain)
  
  SGain = c(0,
          0,1,-1,0,0, # small gain versus neutral
          0,0,0,0,0,0,0,0,0,0
          )
  
  SGain_CueMod = 1/(t(SGain)%*%solve(t(CueMod_mat)%*%as.matrix(CueMod_mat))%*%SGain)
  SGain_AntMod = 1/(t(SGain)%*%solve(t(AntMod_mat)%*%as.matrix(AntMod_mat))%*%SGain)
  SGain_FixMod = 1/(t(SGain)%*%solve(t(FixMod_mat)%*%as.matrix(FixMod_mat))%*%SGain)
  
  
  LGain_bl = c(0,
          1,0,0,0,0, # large gain versus baseline
          0,0,0,0,0,0,0,0,0,0
          )
  
  LGain_bl_CueMod = 1/(t(LGain_bl)%*%solve(t(CueMod_mat)%*%as.matrix(CueMod_mat))%*%LGain_bl)
  LGain_bl_AntMod = 1/(t(LGain_bl)%*%solve(t(AntMod_mat)%*%as.matrix(AntMod_mat))%*%LGain_bl)
  LGain_bl_FixMod = 1/(t(LGain_bl)%*%solve(t(FixMod_mat)%*%as.matrix(FixMod_mat))%*%LGain_bl)
  
  SGain_bl = c(0,
          0,1,0,0,0, # small gain versus baseline
          0,0,0,0,0,0,0,0,0,0
          )
  
  SGain_bl_CueMod = 1/(t(SGain_bl)%*%solve(t(CueMod_mat)%*%as.matrix(CueMod_mat))%*%SGain_bl)
  SGain_bl_AntMod = 1/(t(SGain_bl)%*%solve(t(AntMod_mat)%*%as.matrix(AntMod_mat))%*%SGain_bl)
  SGain_bl_FixMod = 1/(t(SGain_bl)%*%solve(t(FixMod_mat)%*%as.matrix(FixMod_mat))%*%SGain_bl)
  
  
  # VIF w/o intercept
  #diag(solve(cor(CueMod_mat[,2:16])))
  #diag(solve(cor(FixMod_mat[,2:16])))
  #diag(solve(cor(AntMod_mat[,2:16])))
  sub_details = c(sub, run, 
                  LGain_CueMod, LGain_AntMod, LGain_FixMod,
                  SGain_CueMod, SGain_AntMod, SGain_FixMod,
                  # implicit baseline
                  LGain_bl_CueMod, LGain_bl_AntMod, LGain_bl_FixMod,
                  SGain_bl_CueMod, SGain_bl_AntMod, SGain_bl_FixMod)
  
  df = rbind(df,sub_details)
  
}

3.3 Plotting Timings/Efficiency

Pltting the BOLD timeseries based on each cue type from the provided behavioral data. By visualizing the different models, can observe how duration and amplitude shift between them. As a reminder, the models include:

design_CueMod: Cue Onset + Cue Duration
design_AntMod: Cue Onset + (Cue Duration + Fixation Duration)
design_FixMod: Fixation Onset + Fixation Duration

3.3.1 Timing by condition & design

pal <- c("#000000","#004949","#009292","#ff6db6","#ffb6db",
 "#490092","#006ddb","#b66dff","#6db6ff","#b6dbff",
 "#920000","#924900","#db6d00","#24ff24","#ffff6d") 
 
cuemod_plt = design_CueMod %>% 
  gather("Cue","Value", BW:SLmss) %>% 
  ggplot(aes(x = TR, y = Value, colour = Cue)) +
  geom_line() +
  facet_wrap(~Cue) +
  ggtitle("A") +
  #scale_color_brewer(palette = 'Dark2') +
  scale_color_manual(values = pal) +
  theme_minimal()+
  theme(text = element_text(family = "Times New Roman"),
        axis.text = element_text(size = 10),
        axis.title = element_text(size = 10),
        legend.text = element_text(size = 10),
        legend.title = element_text(size = 10),
        plot.title = element_text(size = 16))

antmod_plt = design_AntMod %>% 
  gather("Cue","Value", BW:SLmss) %>% 
  ggplot(aes(x = TR, y = Value, colour = Cue)) +
  geom_line() +
  facet_wrap(~Cue) +
  ggtitle("B") +
  #scale_color_brewer(palette = 'Dark2') +
  scale_color_manual(values = pal) +
  theme_minimal()+
  theme(text = element_text(family = "Times New Roman"),
        axis.text = element_text(size = 10),
        axis.title = element_text(size = 10),
        legend.text = element_text(size = 10),
        legend.title = element_text(size = 10),
        plot.title = element_text(size = 16))
  
fixmod_plt = design_FixMod %>% 
  gather("Cue","Value", BW:SLmss) %>% 
  ggplot(aes(x = TR, y = Value, colour = Cue)) +
  geom_line() +
  facet_wrap(~Cue) +
  ggtitle("C") +
  #scale_color_brewer(palette = 'Dark2') +
  scale_color_manual(values = pal) +
  theme_minimal()+
  theme(text = element_text(family = "Times New Roman"),
        axis.text = element_text(size = 10),
        axis.title = element_text(size = 10),
        legend.text = element_text(size = 10),
        legend.title = element_text(size = 10),
        plot.title = element_text(size = 16))

cuemod_plt

antmod_plt

fixmod_plt

3.3.2 Timing by Design

The below is harder to visualize because the fast event design is shorter than the BOLD signal. Hence, events have quite a bit of overlap.

pal <- c("#000000","#004949","#009292","#ff6db6","#ffb6db",
 "#490092","#006ddb","#b66dff","#6db6ff","#b6dbff",
 "#920000","#924900","#db6d00","#24ff24","#ffff6d") 
 
design_CueMod %>% 
  gather("Cue","Value", BW:SLmss) %>% 
  ggplot(aes(x = TR, y = Value, colour = Cue)) +
  geom_line() +
  ggtitle("Cue Onset + Cue Dur") +
  scale_color_manual(values = pal) +
  theme_minimal() +
  theme(text = element_text(family = "Times New Roman"),
        axis.text = element_text(size = 10),
        axis.title = element_text(size = 10),
        legend.text = element_text(size = 10),
        legend.title = element_text(size = 10),
        plot.title = element_text(size = 16))

design_AntMod %>% 
  gather("Cue","Value", BW:SLmss) %>% 
  ggplot(aes(x = TR, y = Value, colour = Cue)) +
  geom_line() +
  ggtitle("Cue Onset + (Cue Dur + Fix Dur)") +
  scale_color_manual(values = pal) +
  theme_minimal() +
  theme(text = element_text(family = "Times New Roman"),
        axis.text = element_text(size = 10),
        axis.title = element_text(size = 10),
        legend.text = element_text(size = 10),
        legend.title = element_text(size = 10),
        plot.title = element_text(size = 16))

design_FixMod %>% 
  gather("Cue","Value", BW:SLmss) %>% 
  ggplot(aes(x = TR, y = Value, colour = Cue)) +
  geom_line() +
  ggtitle("Fix Onset + Fix Dur") +
  scale_color_manual(values = pal) +
  theme_minimal() +
  theme(text = element_text(family = "Times New Roman"),
        axis.text = element_text(size = 10),
        axis.title = element_text(size = 10),
        legend.text = element_text(size = 10),
        legend.title = element_text(size = 10),
        plot.title = element_text(size = 16))

3.4 Plotting effieincy

Plotting each effieincy model by contrast to show the estimated design efficiency for each datapoint.

# create dataframe of details w/ long format
eff_data = as.data.frame(df)
colnames(eff_data) <- c("Sub","Run",
                        "LarGain_CueMod", "LarGain_AntMod", "LarGain_FixMod",
                        "SmGain_CueMod", "SmGain_AntMod", "SmGain_FixMod",
                        "LarGainvBase_CueMod", "LarGainvBase_AntMod", "LarGainvBase_FixMod",
                        "SmGainvBase_CueMod", "SmGainvBase_AntMod", "SmGainvBase_FixMod")

eff_data_long = eff_data %>% 
  gather("Model","Eff",LarGain_CueMod:SmGainvBase_FixMod) %>%  
  separate(col = Model, into = c("Con","Model"), sep = '_')

eff_data_long$Eff <- as.integer(eff_data_long$Eff)

# plotting
# plot by run and contrast
plot_runcon = eff_data_long %>% ggplot(aes(x = Con, y = Eff, colour = Con, fill = Con)) +
  # setting up raincloud
  geom_rain(alpha = .3, rain.side = 'l',
            boxplot.args.pos = list(
              width = 0.05, position = position_nudge(x = 0.15)),
            violin.args.pos = list(
              side = "l",
              width = 0.7, position = position_nudge(x = 0.3))) +
  # significance comparison
  ggsignif::geom_signif(
  comparisons = list(c("LarGain", "LarGainvBase"),
                     c("LarGainvBase", "SmGain"),
                     c("SmGain","LarGain"),
                     c("SmGainvBase","LarGain")
                     ),
  map_signif_level = TRUE,
  y_position = c(34,31,28,25,22),
  margin_top = 0.1) +
  # plot elements
  theme_classic() +
  scale_fill_brewer(palette = 'Dark2') +
  scale_color_brewer(palette = 'Dark2') +
  guides(fill = 'none', color = 'none') +
  labs(title = "Estimated Efficiency by Run & Contrast",
  caption = "LGain: Large Gain > Neut; SGain: Small Gain > Neut;
  LGain v BL: Large Gain > Implicit Baseline; SGain v BL: Small Gain > Implicit Baseline")+
  ylab("Estimated Effieincy")+
  scale_x_discrete(name = "Model Type",
                   labels = c("LGain","SGain","LGain v BL","SGain v BL"))+
  facet_wrap(~Run + Model)+
  theme(text = element_text(family = "Times New Roman"),
        axis.text = element_text(size = 10),
        axis.title = element_text(size = 10),
        legend.text = element_text(size = 10),
        legend.title = element_text(size = 10),
        plot.title = element_text(size = 16))

plot_runcon

Reliability_SimSummary

¶ Reliability, fMRI, tasks ¶

¶ by Michael Demidenko ¶

March 2023